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Abstract. Recent work on eigenvalues and eigenvectors for tensors of order m > 3 has been 
motivated by applications in blind source separation, magnetic resonance imaging, molecular con- 
formation, and more. In this paper, we consider methods for computing real symmetric-tensor 
eigenpairs of the form 7lx m_1 = Ax subject to ||x|| = 1, which is closely related to optimal rank-1 
__ approximation of a symmetric tensor. Our contribution is a shifted symmetric higher-order power 

,_^ method (SS-HOPM), which we show is guaranteed to converge to a tensor eigenpair. SS-HOPM can 

^"~n be viewed as a generalization of the power iteration method for matrices or of the symmetric higher- 

«i order power method. Additionally, using fixed point analysis, we can characterize exactly which 

eigenpairs can and cannot be found by the method. Numerical examples are presented, including 
^O examples from an extension of the method to finding complex eigenpairs. 
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1. Introduction. Tensor eigenvalues and eigenvectors have received much at- 
tention lately in the literature [13, 18, 20, 19, 4, 15, 27]. The tensor eigenproblem is 
important because it has applications in blind source separation [11], magnetic res- 
onance imaging [25, 22], molecular conformation [7], etc. There is more than one 
possible definition for a tensor eigenpair [18]; in this paper, we specifically use the 
following definition. 

Definition 1.1. Assume that A. is a symmetric m th -order n- dimensional real- 
valued tensor. For any n-dimensional vector x, define 

. n n 

^ (Ax" 1 ^ 1 )^-^--- ^ a ili2 ... im x l2 ---x lm for i 1 = l,...,n. (1.1) 

,— H Then A £1 is an eigenvalue of A if there exists x € M. n such that 

Ax" 1 - 1 = Ax and x T x = 1. (1.2) 



O 



X 



The vector x is a corresponding eigenvector, and (A,x) is called an eigenpair. 

Definition 1.1 is equivalent to the Z-eigenpairs defined by Qi [18, 19] and the 
Z 2 -eigenpairs defined by Lim [13]. In particular, Lim [13] observes that any eigenpair 
(A,x) is a Karush-Kuhn- Tucker (KKT) point (i.e., a constrained stationary point) of 
the nonlinear optimization problem 



Ax m subject to x T x = 1, where ilx m = V' 



max^ix" suDiect to x*x=i, wnere Jtx" = > ••• > a; 

ii— 1 i rn — l 



(1.3) 
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This is equivalent to the problem of finding the best symmetric rank-1 approximation 
of a symmetric tensor [6]. We present the more general definition that incorporates 
complex- valued eigenpairs in §5. 

In this paper, we build upon foundational work by Kofidis and Regalia [11] for 
solving (1.3). Their paper is extremely important for computing tensor eigenvalues 
even though it predates the definition of the eigenvalue problem by three years. Kofidis 
and Regalia consider the higher-order power method (HOPM) [6] , a well-known tech- 
nique for approximation of higher-order tensors, and show that its symmetric gen- 
eralization (S-HOPM) is not guaranteed to converge. They go on, however, to use 
convexity theory to provide theoretical results (as well as practical examples) explain- 
ing conditions under which the method is convergent for even-order tensors (i.e., m 
even). Further, these conditions are shown to hold for many problems of practical 
interest. 

In the context of independent component analysis (ICA) , both Regalia and Kofidis 
[23] and Erdogen [8] have developed shifted variants of the power method and shown 
that they are monotonically convergent. We present a similar method in the context of 
finding real- valued tensor eigenpairs, called the shifted symmetric higher-order power 
method (SS-HOPM), along with theory showing that it is guaranteed to converge to a 
constrained stationary point of (1.3). The proof is general and works for both odd- and 
even-order tensors (i.e., all m > 3). The effectiveness of SS-HOPM is demonstrated 
on several examples, including a problem noted previously [11] for which S-HOPM 
does not converge. We also present a version of SS-HOPM for finding complex-valued 
tensor eigenpairs and provide examples of its effectiveness. 

As mentioned, there is more than one definition of a tensor eigenpair. In the 
case of the l m -eigenpair (we use m for the tensor order instead of k as in some 
references) or H-eigenpair, the eigenvalue equation becomes Ax" 1 ^ 1 = Xxy 71 * 1 ', where 
x [m,— l] denotes the vector x with each element raised to the (to — l) st power [13, 18]. 
In this context, Qi, Wang, and Wang [21] propose some methods specific to third- 
order tensors (to = 3). Unlike the (P-)eigenvalues we consider here, it is possible to 
guarantee convergence to the largest /"-eigenvalue for certain classes of nonnegative 
tensors. For example, see the power methods proposed by Ng, Qi, and Zhou [15] 
and Liu, Zhou, and Ibrahim [14], the latter of which also uses a shift to guarantee 
convergence for any irreducible nonnegative tensor. 

2. Preliminaries. Throughout, let T and £ denote the unit ball and sphere on 
R n , i.e., 

r = {x e R n : ||x|| < 1} and £ = {x e R n : ||xj| = 1}. 

Additionally, define 

n ? „ = the set of all permutations of (1, ... , m). 

Let x_Ly denote x T y = 0, and define x 1 - = {y € R™ : x_Ly}. Let p(A) denote the 
spectral radius of a square matrix A, i.e., the maximum of the magnitudes of its 
eigenvalues. 

2.1. Tensors. A tensor is an m-way array. We let Rl m '™] denote the space of 
TO th -order real-valued tensors with dimension n, e.g., IRl 3,2 ! = IR 2x2x2 . We adopt the 
convention that R [0 ' n] = R. 
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We formally introduce the notion of a symmetric tensor, sometimes also called 
supersymmetric, which is invariant under any permutation of its indices. Further, we 
define a generalization of the tensor-vector multiplication in equations (1.1) and (1.3). 

Definition 2.1 (Symmetric tensor [5]). A tensor A £ R[ m > n l is symmetric if 

a ii-i m for all ii,...,i m £ {1, ...,n} and p £ II m . 



Definition 2.2 (Symmetric tensor-vector multiply). Let A £ R[ m > n l be symmet- 
ric and x £ R". Then for < r < to — I, the (m — r)-times product of the tensor A 
with the vector x is denoted by Ax m_r £ W r,n ' and defined by 

(Ax m ~ r ) il ... ir = ^2 a tl ... im x ir+1 ■■ -x im for all ii, . .. , i r £ {1, . . . ,n}. 

i r+ i,...,i m 

Example 2.3. The identity matrix plays an important role in matrix analysis. 
This notion can be extended in a sense to the domain of tensors. We may define an 
identity tensor as a symmetric tensor £ £ R[ m > n l such that 

Ex" 1 " 1 = x for all x £ S. 

We restrict x £ £ since it is not possible to have a tensor with to > 2 such that the 
above equation holds for all x £ M. n . For any x ^ S, the above equation implies 

Ex™- 1 = ||x|r- 1 £(x/||x||) m - 1 = ||x|| I "- 1 (x/||x||) = ||xH m - 2 x. 

Consider the case of to = 4 and n = 2. The system of equations that must be satisfied 
for all x £ S is 

eimx? + 3em 2 a;?a;2 + 3eii 22 a;ia;2 + ei2 22 a;2 = x u 
eni2^i + 3en 22 a; 1 a; 2 + ie\22iX\x\ + e 2222 £2 = x 2 . 

Consider x = [l 0] . This yields e m i = 1 and em 2 = 0. Similarly, x = [0 l] 
yields e 2222 = 1 and ei 222 = 0. The only remaining unknown is en 22 , and choosing, 
e.g., x = [v2/2 \/2/2] yields en 22 = 1/3. In summary, the identity tensor for 
m = A and n = 2 is 

1 if i = j = k = I, 

1/3 if i = j ^ k = I, 

< ,;i.i = { 1/3 if i = k ^ j = I, 

1/3 iii = l^ j = k, 

otherwise. 

We generalize this idea in the next property. □ 

Property 2.4. For to even, the identity tensor £ £ IR[ m > n ] satisfying £x' Tl_1 = x 
for all x £ S is given by 

6ll '" im ~ ml. 2-~l V(DM2) "*p(3)*p(4) "■■ %(m-l)',(m) (^• 1 ) 

' pen m 
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for ix, . . . ,i m G {1, . . . , n}, where 5 is the standard Kronecker delta, i.e., 

3 \0 ifi^j. 

This identity tensor appears in a previous work [18], where it is denoted by Ie 
and used to define a generalization of the characteristic polynomial for symmetric 
even-order tensors. 

Example 2.5. There is no identity tensor for m odd. This is seen because 
if £x m_1 = x for some odd m and some x G E, then we would have x G E but 
£(-x) m - 1 =x^ -x. D 

For any even-order tensor (i.e., to even), observe that if (A, x) is an eigenpair, 
then (A, — x) is also an eigenpair since 

A(-x) m - 1 = -Ax" 1 " 1 = A(-x). 

Likewise, for any odd-order tensor (i.e., m odd), (—A, — x) is also an eigenpair since 

A(-x)™- 1 = Ax™" 1 = (-A)(-x). 

These are not considered to be distinct eigenpairs. 

We later present, as Theorem 5.3, a recently derived result [3] that bounds the 
number of real eigenpairs by ((to— 1)™ — l)/(m— 2). We defer discussion of this result 
until §5, where we discuss complex eigenpairs. 

Because the tensor eigenvalue equation for m > 2 amounts to a system of non- 
linear equations in the components of x, a direct solution is challenging. Numerical 
algorithms exist for finding all solutions of a system of polynomial equations, but be- 
come computationally expensive for systems with many variables (here, large n) and 
with high-order polynomials (here, large to). A polynomial system solver (NSolve) 
using a Grobner basis method is available in Mathematica [28] and has been employed 
to generate a complete list of eigenpairs for some of the examples in this paper. The 
solver is instructed to find all solutions (A, x) of the system (1.2). Redundant solutions 
with the opposite sign of x (for even to) or the opposite signs of x and A (for odd to) 
are then eliminated. 

2.2. Convex functions. Convexity theory plays an important role in our anal- 
ysis. Here we recall two important properties of convex functions [2] . 

Property 2.6 (Gradient of convex function). A differentiable function f : fl C 
W 1 — > R is convex if and only if ft is a convex set and /(y) > /(x) + V/(x) T (y — x) 
for all x, y G Q. 

Property 2.7 (Hessian of convex function). A twice differentiable function 
/:1)C R™ — ^ R is convex if and only if fl is a convex set and the Hessian 1 of f is 
positive semidefinite on fi, i.e., V 2 /(x) y for all x e Q. 

We prove an interesting fact about convex functions on vectors of unit norm that 
will prove useful in our later analysis. This fact is implicit in a proof given previously 
[11, Theorem 4] and explicit in [23, Theorem 1]. 

Theorem 2.8 (Kofidis and Regalia [11, 23]). Let f be a function that is con- 
vex and continuously differentiable on T. Let w G E with V/(w) 7^ 0. If v = 
V/(w)/||V/(w)|| ^ w, then /(v) - /(w) > 0. 



x By V 2 we denote the Hessian matrix and not its trace, the Laplacian. 
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Proof. For arbitrary nonzero z G W l , z x is strictly maximized for x £ E by 
x = z/||z||. Substituting z = V/(w), it follows that V/(w) T v > V/(w) T w, since 
v = V/(w)/||V/(w)|| 7^ w and w G E. By the convexity of / on T and Property 2.6, 
we have /(v) > /(w) + V/(w) T (v — w) for all v, w G T. Consequently, /(v) — /(w) > 
V/(w) T (v-w) >0. □ 

2.3. Constrained optimization. Here we extract relevant theory from con- 
strained optimization [17]. 

Theorem 2.9. Let f : M™ — > E 6e continuously differentiable. A point x* G E is 
a (constrained) stationary point of 

max/(x) subject to x G E 

i/ t/iere exists /i* G E suc/i t/iat V/(x„) + /i*x* = 0. TTie point x* is a (constrained) 
isolated local maximum if additionally, 

w T (V 2 /(x,)+^J)w < /or s/i weEHx^. 

Proof. The constraint x G E can be expressed as c(x) = ^(x T x — 1) = 0. The 
Lagrangian for the constrained problem is then given by 

£(x,M) = /(*) +M C ( X )- 

Its first and second derivatives with respect to x are 

V£(x, (i) = V/(x) + M x and V 2 £(x, ^) = V 2 /(x) + /il. 

By assumption, V£(x*,/i*) = and c(x*) = 0. Therefore, the pair (x*,//*) satisfies 
the Karush-Kuhn- Tucker (KKT) conditions [17, Theorem 12.1] and so is a constrained 
stationary point. It is additionally a constrained isolated local maximum if it meets 
the second-order sufficient condition [17, Theorem 12.6]. D 

2.4. Fixed point theory. We consider the properties of iterations of the form 

Xfc+l = 0(x fe ). 

Under certain conditions, the iterates are guaranteed to converge to a fixed point. In 
particular, we are interested in "attracting" fixed points. 

Definition 2.10 (Fixed point). A point x* G W 1 is a fixed point of cf> : W l -> R™ 
if 4>(x.*) = x*. Further, x* is an attracting fixed point if there exists 5 > such that 
the sequence {xj,} defined by Xfc +1 = 0(xj.) converges to x* for any x such that 
||x -x*|| < S. 

Theorem 2.11 ([24, Theorem 2.8]). Letx* G E™ be afixedpoint of<f> : E™ -> E n ; 
and Zet J : R n — >• R" xn oe i/ie Jacobian of (f). Then x* is an attracting fixed point if 
a = p(J(x*)) < 1; further, if a > 0, then the convergence of Xk+i — </>(xfc) to x* is 
linear with rate a. 

This condition on the Jacobian for an attracting fixed point is sufficient but not 
necessary. In particular, if a = /j(«7(x*)) = 1, then x* may or may not be attracting, 
but there is no neighborhood of linear convergence to it. For a < 1, the rate of linear 
convergence depends on a and is slower for a values closer to 1. On the other hand, 
for a > 1, an attractor is ruled out by the following. 

Theorem 2.12 ([26, Theorem 1.3.7]). Let x„ G E" be a fixed point of <f> : E" -> 
E n , and let J : R n — > R raxn be the Jacobian of <f>. Then x* is an unstable fixed point 
if o~ = yo(J(x*)) > 1. 
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3. Symmetric higher-order power method (S-HOPM). We review the 
symmetric higher-order power method (S-HOPM), introduced by De Lathauwer et 
al. [6] and analyzed further by Kofidis and Regalia [11]. The purpose of S-HOPM is 
to solve the optimization problem 

max|.Ax m | subject to x G £. (3.1) 

The solution of this problem will be a solution of either the following maximization 
problem (lacking the absolute value) or its opposite minimization problem: 

max /(x) subject to x G S, where /(x) = Ax m . (3-2) 

Setting A = /(x), these problems are equivalent to finding the best symmetric rank-1 
approximation of a symmetric tensor A G lS} m ' n \ i.e., 

min \\A — S|| subject to 6i 1 ...i m = Ax^ • • • xi m and x G S. (3-3) 

A.x 

Details of the connection between (3.2) and (3.3) are available elsewhere [6]. The S- 
HOPM algorithm is shown in Algorithm 1. We discuss its connection to the eigenvalue 
problem in §3.1 and its convergence properties in §3.2. 

Algorithm 1 Symmetric higher-order power method (S-HOPM) [6, 11] 
Given a symmetric tensor A G ]R[ m '"l. 
Require: x G E n with ||x || = 1. Let A = Axff. 
l: for k = 0, 1, ... do 



x fc+1 <- .Ax™- 1 

X fc+ i <r- X fc+ l/||x fc+ i| 

A/c+i «- Ax™ +1 
end for 



3.1. Properties of /(x) = Ax." 1 . The function f(x) = Ax™ 1 plays an important 
role in the analysis of eigenpairs of A because all eigenpairs are constrained stationary 
points of /, as we show below. 

We first need to derive the gradient of /. This result is perhaps generally well 
known [13, Equation 4], but here we provide a proof. 

Lemma 3.1. Let A G Rl m ' n l be symmetric. The gradient o//(x) = Ax m is 

9 (x)EV/(x)=mAx m - 1 eR". (3.4) 

Proof. We use the basic relation S7k%j = Sjk- Applying the product rule to (3.2), 
we find 

m 
il,...,i m q— 1 

Upon bringing the sum over q to the outside, we observe that for each q the dummy 
indices i± and i q can be interchanged (without affecting the symmetric tensor A), and 
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the result is independent of q: 

m 
Vfc/^Xj = > > CLi 1 i 2 ...i rn Oi 1 kXi 2 ■ • • Xi q _ 1 Xi q Xi q + 1 • • • Xi m 

q—1 ii,...,i m 
m 

q—1 i 2 , ■■■,im 

= m(Ax™- 1 ),.. 

Hence, V/(x) = rnXx"*- 1 . D 

Theorem 3.2. Let A € Rl m >"J be symmetric. Then (A, x) is an eigenpair of A 
if and only if x is a constrained stationary point of (3.2). 

Proof. By Theorem 2.9, any constrained stationary point x* of (3.2) must satisfy 
mTlx" 1-1 + jU*x* = for some /i* £ K. Thus, A* = — /z*/m is the eigenvalue corre- 
sponding to x* . Conversely, any eigenpair meets the condition for being a constrained 
stationary point with fj, t = — toA*. □ 

This is is the connection between (3.2) and the eigenvalue problem. It will also 
be useful to consider the Hessian of /, which we present here. 

Lemma 3.3. Let A e Rl m,n l be symmetric. The Hessian o//(x) = ilx"' is 

H(x) = V 2 /(x) = m(m - lJAx" 1 " 2 e R nxn . (3.5) 

Proof. The (j, k) entry of -ff(x) is given by the k th entry of V<?j(x). The function 
<7j(x) can be rewritten as 

g j (x) = m ^2 c >'j i2 ...i m Xi 2 ---Xi m = m^B^'x m ~ 1 



where "B^' is the order-(m — 1) symmetric tensor that is the j th subtensor of A, 

.,■ , . From Lemma 3.1, we he 

V&(x) = m(m - l)2 (j) x m - 2 . 



defined by 6„ - =a„- 1 ...,- , . From Lemma 3.1, we have 



Consequently, 



(H(x))j k = m(m - 1) J^ ajki 3 -i m x is ■ ■ ■ x 



lrn ? 



that is, ff(x) = to(to - l)Ax'"- 2 . D 

From Theorem 2.9, we know that the projected Hessian of the Lagrangian plays 
a role in determining whether or not a fixed point is a local maximum or minimum. 
In our case, since [i* = — mA*, for any eigenpair (A*,x*) (which must correspond to 
a constrained stationary point by Theorem 3.2) we have 

V 2 £(x,, A,) = m(m - I) Ax", 1 - 2 - toAJ. 

Specifically, Theorem 2.9 is concerned with the behavior of the Hessian of the La- 
grangian in the subspace orthogonal to x*. Thus, we define the projected Hessian of 
the Lagrangian as 

C(A*,x*) = U*" ((to - l)Ax™- 2 - A*I) U* e r(™-i)x(«-i) ; (3.6) 



8 T. G. Kolda and J. R. Mayo 

where the columns of U* € jj™*!.™- 1 ) form an orthonormal basis for x^J". Note that 
we have removed a factor of m for convenience. We now classify eigenpairs according 
to the spectrum of C(A*,x*). The import of this classification will be made clear in 
§4.2. 

Definition 3.4. Let A € W- m>n ' be a symmetric tensor. We say an eigenpair 
(A,x) of A € Ml" 1 -™! is positive stable ifC(\, x) is positive definite, negative stable if 
C(A, x) is negative definite, and unstable if C(A, x) is indefinite. 

These labels are not exhaustive because we do not name the cases where C(A, x) is 
only semidefinite, with a zero eigenvalue. Such cases do not occur for generic tensors. 

If m is odd, then (A, x) is positive stable if and only if (—A, — x) is negative stable, 
even though these eigenpairs are in the same equivalence class. On the other hand, 
if m is even, then (A, x) is a positive (negative) stable eigenpair if and only if (A, — x) 
is also positive (negative) stable. 

3.2. S-HOPM convergence analysis. S-HOPM has been deemed unreliable 
[6] because convergence is not guaranteed. Kofidis and Regalia [11] provide an analysis 
explaining that S-HOPM will converge if certain conditions are met, as well as an 
example where the method does not converge, which we reproduce here. 

Example 3.5 (Kofidis and Regalia [11, Example 1]). Let A e K [4 ' 3] be the 
symmetric tensor defined by 

ami = 0.2883, a m2 = -0.0031, a m3 = 0.1973, 01122 = -0.2485, 

01123 = -0.2939, 01133 = 0.3847, a 1222 = 0.2972, a 1223 = 0.1862, 

a 1233 = 0.0919, a 1333 = -0.3619, a 2222 = 0.1241, a 2223 = -0.3420, 

a 2 233- 0.2127, a 2333 = 0.2727, a 3333 = -0.3054. 

Kofidis and Regalia [11] observed that Algorithm 1 does not converge for this tensor. 
Because this problem is small, all eigenpairs can be calculated by Mathematica as 
described in §2.1. From Theorem 5.3, this problem has at most 13 eigenpairs; we 
list the 11 real eigenpairs in Table 3.1. We ran 100 trials of S-HOPM using differ- 
ent random starting points xo chosen from a uniform distribution on [—1, 1]". For 
these experiments, we allow up to 1000 iterations and say that the algorithm has 
converged if |A fe+1 — A fc | < 10~ 16 . In every single trial for this tensor, the algo- 
rithm failed to converge. In Figure 3.1, we show an example {A^} sequence with 
xo = [—0.2695 0.1972 0.3370] . This coincides with the results reported previ- 
ously [11]. □ 

Example 3.6. As a second illustrative example, we consider an odd-order tensor 
A € M [3 ' 3] defined by 



8111 = 


-0.1281, 


GS112 = 


0.0516, 


0113 = 


-0.0954, 


0-122 — 


-0.1958, 


Ol23 = 


-0.1790, 


a 133 = 


-0.2676, 


a 222 = 


0.3251, 


0-223 = 


0.2513, 


«233 = 


0.1773, 


a 333 = 


0.0338. 











From Theorem 5.3, A has at most 7 eigenpairs; in this case we achieve that bound and 
the eigenpairs are listed in Table 3.2. We ran 100 trials of S-HOPM as described for 
Example 3.5. Every trial converged to either A = 0.8730 or A = 0.4306, as summarized 
in Table 3.3. Therefore, S-HOPM finds 2 of the 7 possible eigenvalues. □ 

In their analysis, Kofidis and Regalia [11] proved that the sequence {\k} in Al- 
gorithm 1 converges if A € JRl™ 1 -™] is even-order and the function /(x) is convex or 
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Table 3.1: Eigenpairs for A £ B^ 4,3 ' from Example 3.5. 



A 


x r 


Eigenvalues of C(A, x) 


Type 


0.8893 


[ 0.6672 0.2471 -0.7027 ] 


{ -0.8857, -1.8459 } 


Ncg. stable 


0.8169 


0.8412 -0.2635 0.4722 ] 


{ -0.9024, -2.2580 } 


Neg. stable 


0.5105 


0.3598 -0.7780 0.5150 ] 


{ 0.5940, -2.3398 } 


Unstable 


0.3633 


0.2676 0.6447 0.7160 ] 


{ -1.1765, -0.5713 } 


Neg. stable 


0.2682 


0.6099 0.4362 0.6616 ] 


{ 0.7852, -1.1793 } 


Unstable 


0.2628 


0.1318 -0.4425 -0.8870 ] 


{ 0.6181, -2.1744 } 


Unstable 


0.2433 


0.9895 0.0947 -0.1088 ] 


{ -1.1942, 1.4627 } 


Unstable 


0.1735 


1 0.3357 0.9073 0.2531 ] 


{ -1.0966, 0.8629 } 


Unstable 


-0.0451 


0.7797 0.6135 0.1250 ] 


{ 0.8209, 1.2456 } 


Pos. stable 


-0.5629 


0.1762 -0.1796 0.9678 ] 


{ 1.6287, 2.3822 } 


Pos. stable 


-1.0954 


0.5915 -0.7467 -0.3043 ] 


{ 1.8628, 2.7469 } 


Pos. stable 



^ 




20 40 60 80 100 

Iteration (k) 



Fig. 3.1: Example A fc values for S-HOPM on A £ Rl 4 ' 3 ! from Example 3.5. 



Table 3.2: Eigenpairs for A £ 



>[3,3] 



from Example 3.6. 



A 


x J 


Eigenvalues of C(A, x) 


Type 


0.8730 


[ -0.3922 0.7249 0.5664 ] 


{ -1.1293, -0.8807 } 


Neg. stable 


0.4306 


-0.7187 -0.1245 -0.6840 ] 


{ -0.4420, -0.8275 } 


Neg. stable 


0.2294 


-0.8446 0.4386 -0.3070 ] 


{ -0.2641, 0.7151 } 


Unstable 


0.0180 


0.7132 0.5093 -0.4817 ] 


{ -0.4021, -0.1320 } 


Neg. stable 


0.0033 


0.4477 0.7740 -0.4478 ] 


{ -0.1011, 0.2461 } 


Unstable 


0.0018 


0.3305 0.6314 -0.7015 ] 


{ 0.1592, -0.1241 } 


Unstable 


0.0006 


0.2907 0.7359 -0.6115 ] 


{ 0.1405, 0.0968 } 


Pos. stable 



Table 3.3: Eigenpairs for A £ 
100 random starts. 



»[3,3] f rom Example 3.6 computed by S-HOPM with 



# Occurrences 


A 


X 


Median Its. 


62 


0.8730 


[ -0.3922 0.7249 0.5664 ] 


19 


38 


0.4306 


-0.7187 -0.1245 -0.6840 ] 


184 
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concave on M. n . Since m — 2£ (because m is even), / can be expressed as 
/(x) = (x(g)-- -(8ix ) T A( x(8i ■ ■ •®x ), 

& times £ times 

where A £ M. n xn is an unfolded version of the tensor A. 2 Since A is symmetric, it 
follows that A is symmetric. The condition that / is convex (concave) is satisfied if 
the Hessian 

V 2 /(x) = (I® x® ■ ■ •®x ) T A(I® x® ■ • •®x ) 

I — 1 times £ — 1 times 

is positive (negative) semidefinite for all x £ W 1 . 

We make a few notes regarding these results. First, even though / is convex, its 
restriction to the nonconvex set E is not. Second, {Afc} is increasing if / is convex 
and decreasing if / is concave. Third, only {A&} is proved to converge for S-HOPM 
[11, Theorem 4]; the iterates {x^} may not. In particular, it is easy to observe that 
the sign of x^ may flip back and forth if the concave case is not handled correctly. 

4. Shifted symmetric higher-order power method (SS-HOPM). In this 
section, we show that S-HOPM can be modified by adding a "shift" that guarantees 
that the method will always converge to an eigenpair. In the context of ICA, this 
idea has also been proposed by Regalia and Kofidis [23] and Erdogen [8]. Based on 
the observation that S-HOPM is guaranteed to converge if the underlying function is 
convex or concave on K n , our method works with a suitably modified function 

/(x) = /(x) + a(x T x)"V 2 . (4.1) 

Maximizing / on E is the same as maximizing / plus a constant, yet the properties 
of the modified function force convexity or concavity and consequently guarantee 
convergence to a KKT point (not necessary the global maximum or minimum). Note 
that previous papers [23, 8] have proposed similar shifted functions that are essentially 
of the form /(x) = /(x) + ax T x, differing only in the exponent. 

An advantage of our choice of / in (4.1) is that, for even m, it can be interpreted 



/(x) = Ax m = (A + a£)x m , 

where £ is the identity tensor as defined in (2.1). Thus, for even m, our proposed 
method can be interpreted as S-HOPM applied to a modified tensor that directly 
satisfies the convexity properties to guarantee convergence [11]. Because £x m_1 = x 
for x £ E, the eigenvectors of A are the same as those of A and the eigenvalues are 
shifted by a. Our results, however, are for both odd- and even-order tensors. 

Algorithm 2 presents the shifted symmetric higher-order power method (SS- 
HOPM). Without loss of generality, we assume that a positive shift (a > 0) is used 
to make the modified function in (4.1) convex and a negative shift (a < 0) to make 
it concave. We have two key results. Theorem 4.4 shows that for any starting point 



2 Specifically, A = A-iizxC) with 1Z = {1. . . . ,£} and C = {£ + 1, . . . , m} in matricization notation 

[12]. 
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xo £ E, the sequence {A^} produced by Algorithm 2 is guaranteed to converge to an 
eigenvalue in the convex case if 

a > /3(A) = (m - 1) ■ maxp(Ax m " 2 ). (4.2) 

x£E 

Corollary 4.6 handles the concave case where we require a < —/3(A). Theorem 4.8 
further shows that Algorithm 2 in the convex case will generically converge to an 
eigenpair (A, x) that is negative stable. Corollary 4.9 proves that Algorithm 2 in the 
concave case will generically converge to an eigenpair that is positive stable. Generally, 
neither version will converge to an eigenpair that is unstable. 

Algorithm 2 Shifted Symmetric Higher-Order Power Method (SS-HOPM) 
Given a tensor A £ E. [m - n] . 

Require: x £ M. n with ||x || = 1. Let A = Axff. 
Require: a£R 

l: for k = 0, 1, . . . do 

2: if a > then 

3: Xfe + i <— ybc™ -1 + axfe > Assumed Convex 

4: else 

5: Xfe + i i (Ax™ -1 + axfc) > Assumed Concave 

6: end if 

7: X fc+1 <- X/. + i/||x fe+ i|| 

8: A fe+ i «- Ax% +1 

9: end for 



4.1. SS-HOPM convergence analysis. We first establish a few key lemmas 
that guide the choice of the shift a > (3(A) in SS-HOPM. 

Lemma 4.1. Let A £ R[ m >"] be symmetric and let /3(A) be as defined in (4-2). 
Then /3(A) < (to- 1)^,...^ K....J- 

Proof. For all x, y £ E, we obtain |y T (Ax'"" 2 )y| < ^h,...,^ I a ii...i m I b y applying 
the triangle inequality to the sum of n m terms. Thus p(Ax. m ~ 2 ) < J^i i \ a ii...i m \ 
for all x £ E, and the result follows. D 

Lemma 4.2. Let A £ Rl m '"l be symmetric, let /(x) = Ax™, and let [3(A) be as 
defined in (4.2). Then |/(x)| < (3(A) /(m - 1) for all x £ E. 

Proof. We have |.Ax m | = |x T (Ax m ~ 2 )x| < p(Ax m - 2 ) < /3(A)/(m- 1). D 

The preceding lemma upper bounds the magnitude of any eigenvalue of A by 
(3(A) /(m— 1) since any eigenpair (A,x) satisfies A = /(x). Thus, choosing a > (3(A) 
implies that a is greater than the magnitude of any eigenvalue of A. 

Lemma 4.3. Let A £ R[ m > n l be symmetric and let H(x) and (3(A) be as defined 
in (3.5) and (4.2). Then p(H(x)) < m/3(A) for all x £ E. 

Proof. This follows directly from (3.5) and (4.2). D 

The following theorem proves that Algorithm 2 will always converge. Choosing 
a > (to — 1) Si i \ a ii--im\ i s a conservative choice that is guaranteed to work 
by Lemma 4.1, but this may slow down convergence considerably, as we show in 
subsequent analysis and examples. 

Theorem 4.4. Let A £ Rl" 1 -"] be symmetric. For a > (3(A), where (3(A) is 
defined in (4-2), the iterates {Afe,Xfe} produced by Algorithm 2 satisfy the following 
properties, (a) The sequence {Xk} is nondecreasing, and there exists A* such that 
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Xk — > A*, (b) The sequence {xfc} has an accumulation point, (c) For every such 
accumulation point x*, the pair (A*,X*) is an eigenpair of A., (d) If A has finitely 
many real eigenvectors, then there exists x* such that Xfc — > x* . 

Proof. Our analysis depends on the modified function / defined in (4.1). Its 
gradient and Hessian for x^O are 

g(x) = V/(x) = g(x) + mafx^)" 1 ' 2 -^, (4.3) 

H(x) = V 2 f(x) - H(x) + ma(x T x) m / 2 - 1 I + m(m - 2)a(x T x) m ^- 2 xx T , (4.4) 

where g and H are the gradient and Hessian of / from Lemma 3.1 and Lemma 3.3, 
respectively. And because f(x) = 0(||x|j m ), as x — > 0, it follows that /(x) is of third 
or higher order in x for m > 3; thus g(0) — and H(0) = 0. 

Because it is important for the entire proof, we first show that / is convex on R™ 
for a > (3(A). As noted, if x = 0, we have H(x) = for m > 3. Consider nonzero 
x £ M. n and define x = x/||x|| £ E; then -ff(x) is positive semidefinite (in fact, positive 
definite) by Lemma 4.3 since 

y T H(x)y = ||x|| m " 2 (y T H(x)y + ma + m(m - l)a(x T y) 2 ) 

> ||x|| m " 2 (-m/3(A)+mQ + 0) > for all y e S. 

By Property 2.7, / is convex on M. n because its Hessian is positive semidefinite. 

We also note that —a must be an eigenvalue of A if g(x) = for some x € S, 
since 

g(x) = implies Ax" 1 ^ 1 +ax = 0. 

By Lemma 4.2, choosing a > (1(A) ensures that a is greater than the magnitude 
of any eigenvalue, and so g(x) ^ for all x e E. This ensures that the update in 
Algorithm 2, which reduces to 

x fc+1 = j^pXu (4-5) 

in the convex case, is always well defined. 

(a) Since / is convex on L and x fe+ i,x fe G E and x fe+ i = V/(x fc )/||V/(x fc )||, 
Theorem 2.8 yields 

Afe+i - A fe = f(x k+ i) - f(x k ) > 0, 

where the nonstrict inequality covers the possibility that x/j + i = x k . Thus, {Afc} is 
a nondecreasing sequence. By Lemma 4.2, Afc = f(x k ) is bounded, so the sequence 
must converge to a limit point A*. 3 

(b) Since {xfc} is an infinite sequence on a compact set E, it must have an accu- 
mulation point x* G E by the Bolzano- Weierstrass theorem. Note also that continuity 
of / implies that A* = Ax™. 

(c) By part (a) of the proof, convexity of /, and Property 2.6, we have 

/( X fc +1 ) - / (xfc) -> 



3 Note that the similar approach proposed for ICA [23, Theorem 2] allows the shift a to vary at 
each iteration so long as the underlying function remains convex. 
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and thus 

.g(x fe ) T (x fe+1 - x fc ) -> 0. 
Using (4.5), we can rewrite the above formula as 

||g(x fc )||-g(x fe ) T x fe ^0. (4.6) 

By continuity of g, an accumulation point x* must satisfy 

||.g(x >f )||-.g(x,) T x >f =0, (4.7) 

which implies 

||<7(x*)|| = g(x*) T x* = (rnAx" 1 ^ 1 + max*) 7 *.* = m(A* + a). 

Because x* £ X, (4.7) can hold only if 

<?(x„) mAx™^ 1 + max* 
X * = 1x1 = m(K+a) ' 

that is, 

Hence (A*,x*) is an eigenpair of A. 
(d) Equation (4.6) gives 

||5(x fe )||(l-x£ +1 x fc )^0. 

Because ||g(xfe)|| is bounded away from and because x^, x^+i £ S, this requires that 

Hxfc-Xfc+ill-^0. (4.8) 

Recall that every accumulation point of {x^} must be a (real) eigenvector of A. If 
these eigenvectors are finite in number and thus isolated, consider removing an arbi- 
trarily small open neighborhood of each from E, leaving a closed and thus compact 
space Y C S containing no accumulation points of {x^.}. If {x^.} had infinitely many 
iterates in 7, it would have an accumulation point in Y by the Bolzano- Weierstrass 
theorem, creating a contradiction. Therefore at most finitely many iterates are in Y, 
and {xj.} is ultimately confined to arbitrarily small neighborhoods of the eigenvec- 
tors. By (4.8), however, ||x/- — Xfc + i|| eventually remains smaller than the minimum 
distance between any two of these neighborhoods. Consequently, the iteration ulti- 
mately cannot jump from one neighborhood to another, and so in the limit {x^} is 
confined to an arbitrarily small neighborhood of a single eigenvector x* , to which it 
therefore converges. 
Hence, the proof is complete. D 

Note that the condition of finitely many real eigenvectors in part (d) holds for 
generic tensors. We conjecture that the convergence of {x^} is guaranteed even with- 
out this condition. 

Example 4.5. Again consider A £ JRl 4,3 ' from Example 3.5. We show results us- 
ing a shift of a = 2. We ran 100 trials of SS-HOPM using the experimental conditions 
described in Example 3.5. We found 3 real eigenpairs; the results are summarized in 
Table 4.1a. Three example runs (one for each eigenvalue) are shown in Figure 4.1a. 
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We also considered the "conservative" choice ol a = (rn — 1) J2i i l a ii...i m | — 
55.6620. We ran 100 trials of SS-HOPM using the experimental conditions described 
in Example 3.5, except that we increased the maximum number of iterations to 10,000. 
Every trial converged to one of the same 3 real eigenpairs, but the number of iterations 
was around 1000 (versus around 60 for a = 2); in §4.2, we see that the rate of 
convergence asymptotically decreases as a increases. 

Analogous results are shown for A. € R[ 3,3 J from Example 3.6 with a shift of a = 1 
in Table 4.1b and Figure 4.1b. Here SS-HOPM finds 2 additional eigenpairs compared 
to S-HOPM. In this case, we also considered a — (m — 1) J^\ i \cn x ...i m \ = 9.3560, 
but this again increased the number of iterations up to a factor of ten. 

For both tensors, {Xk} is always a nondecreasing sequence. Observe further that 
SS-HOPM converges only to eigenpairs that are negative stable. □ 



# Occurrences 


A 


X 


Median Its. 


46 


0.8893 


[ 0.6672 


0.2471 


-0.7027 ] 


63 


24 


0.8169 


0.8412 


-0.2635 


0.4722 ] 


52 


30 


0.3633 


0.2676 


0.6447 


0.7160 ] 


65 



(a) A 6 M' 4 ' 3 ' from Example 3.5 with . 



# Occurrences 


A 


X 


Median Its. 


40 


0.8730 


[ -0.3922 


0.7249 


0.5664 ] 


32 


29 


0.4306 


[ -0.7187 


-0.1245 


-0.6840 ] 


48 


18 


0.0180 


0.7132 


0.5093 


-0.4817 ] 


116 


13 


-0.0006 


[ -0.2907 


-0.7359 


0.6115 ] 


145 



(b) A G IR [3 ' 31 from Example 3.6 with a = 1. 

Table 4.1: Eigenpairs computed by SS-HOPM (convex) with 100 random starts. 



SS-HOPM with a = 2 



SS-HOPM with < 




-0.5 



40 60 80 100 

Iteration (k) 

(a) A £ M [4,3 l from Example 3.5 with a = 2. 



05 



-05 



20 













r~ 




-0 0006 


L^ 




¥ 




y~~ 



80 



100 



40 60 

Iteration (k) 

(b) A e K [3 ' 3 ' from Example 3.6 with a = 1 



Fig. 4.1: Example A& values for SS-HOPM (convex). One sequence is shown for each 
distinct eigenvalue. 



Using a large enough negative value of a makes / concave. It was observed 
[11] that /(x) = /(— x) for even-order tensors and so the sequence {Xk} converges 



Shifted Power Method for Computing Tensor Eigenpairs 



15 



regardless of correctly handling the minus sign. The only minor problem in the concave 
case is that the sequence of iterates {x^} does not converge. This is easily fixed, 
however, by correctly handling the sign as we do in Algorithm 2. The corresponding 
theory for the concave case is presented in Corollary 4.6. In this case we choose a to 
be negative, i.e., the theory suggests a < —(3(A). 

Corollary 4.6. Let A e Rt m >"] be symmetric. For a < -f3(A), where /3(A) 
is defined in (4-2), the iterates {Afc,Xfc} produced by Algorithm 2 satisfy the following 
properties, (a) The sequence {Xk} is nonincreasing , and there exists A* such that 
Afc — > A*, (b) The sequence {x^,} has an accumulation point, (c) For any such 
accumulation point x* ; the pair (A*,X*) is an eigenpair of A. (d) If the eigenvalues 
of A. are isolated, then x^ — > x* . 

Proof. Apply the proof of Theorem 4.4 with /(x) = —Ax m . D 

Example 4.7. Revisiting A e Ml 4,3 ' in Example 3.5 again, we run another 100 
trials using a = —2. We find 3 (new) real eigenpairs; the results are summarized in 
Table 4.2a. Three example runs (one for each eigenvalue) are shown in Figure 4.2a. 

We also revisit A £ Ml 3,3 ! from Example 3.6 and use a = — 1. In this case, we 
find the opposites, i.e., (— A,— x), of the eigenpairs found with a = 1, as shown in 
Table 4.2b. This is to be expected for odd-order tensors since there is symmetry, i.e., 
/(x) = — /(— x), C(A,x) = — C{— A,— x), etc. Observe that the median number of 
iterations is nearly unchanged; this is explained in the subsequent subsection where 
we discuss the rate of convergence. Four example runs (one per eigenvalue) are shown 
in Figure 4.2b. 

The sequence {Xk} is nonincreasing in every case. Each of the eigenpairs found 
in the concave case is positive stable. □ 



# Occurrences 


A 


X 


Median Its. 


15 


-0.0451 


[ -0.7797 -0.6135 -0.1250 ] 


35 


40 


-0.5629 


[ -0.1762 0.1796 -0.9678 ] 


23 


45 


-1.0954 


-0.5915 0.7467 0.3043 ] 


23 



(a) A G R [4 ' 3] from Example 3.5 with a = -2. 



# Occurrences 


A 


X 


Median Its. 


19 


0.0006 


[ 0.2907 0.7359 -0.6115 ] 


146 


18 


-0.0180 


-0.7132 -0.5093 0.4817 ] 


117 


29 


-0.4306 


0.7187 0.1245 0.6840 ] 


49 


34 


-0.8730 


0.3922 -0.7249 -0.5664 ] 


33 



(b) A e Ml 3 ' 3 ' from Example 3.6 with a = -1. 

Table 4.2: Eigenpairs computed by SS-HOPM (concave) with 100 random starts. 



4.2. SS-HOPM fixed point analysis. In this section, we show that fixed 
point analysis allows us to easily characterize convergence to eigenpairs according to 
whether they are positive stable, negative stable, or unstable. The convex version of 
SS-HOPM will generically converge to eigenpairs that are negative stable; the concave 
version will generically converge to eigenpairs that are positive stable. 

To justify these conclusions, we consider Algorithm 2 in the convex case as a fixed 
point iteration Xfc + i = 0(x^; a), where <fi is defined as 



(x; a) = 0i(</> 2 (x; a)) with 0i(x) 



(x T x)t 



and (f>2 (x; a) = Ax 1 ' 



(4.9) 
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SS-HOPM with a = -2 



SS-HOPM with a = -1 



0.5 









n rmc-i 








I 

I 


- 



0.5 



-0.5 



0.0006 
-0.0180 
-0.4306 
-0.8730 



^~ 




20 



80 



100 



40 60 

Iteration (k) 

(a) A £ Ml 4 ' 3 ) from Example 3.5 with a = 2. 



20 



80 



100 



40 60 

Iteration (k) 

(b) A e Rl 3 ' 3 ' from Example 3.6 with a = 1 



Fig. 4.2: Example Afc values for SS-HOPM (concave). One sequence is shown for each 
distinct eigenvalue. 



Note that an eigenpair (A, x) is a fixed point if and only if A + a > 0, which is always 
true for a > f3(A). 

From [9] , the Jacobian of the operator <fi is 

J(x;a) = 0iO 2 (x;a))02(x;a), 

where derivatives are taken with respect to x and 

#(x)= (XTX) ^ 3 XXT and «/> 2 (x;a) = (m-l)Ax m - 2 + aI. 
(x-* x) 2 

At any eigenpair (A,x), we have 



2 (x;a) = (A + a)x, ^i(^ 2 (x;a)) 



A + a 
and 2 (x; a) = (m - l)Ax m_2 + al. 



Thus, the Jacobian at x is 



J(x; a) = 



(m - l)(Ax m - 2 - Axx J ) + a(I - xx J ) 
X + a 



(4.10) 



Observe that the Jacobian is symmetric. 

Theorem 4.8. Let (A, x) be an eigenpair of a symmetric tensor A. £ My n ' n '. 
Assume a £ M. such that a > /8(A), where /3(A) is as defined in (4-2). Let </>(x) be 
given by (4-9). Then (A,x) is negative stable if and only if x is a linearly attracting 
fixed point of cj). 

Proof. Assume that (A, x) is negative stable. The Jacobian J(x; a) is given by 
(4.10). By Theorem 2.11, we need to show that p(J(x;a)) < 1 or, equivalently since 
Jfaa) is symmetric, |y T J(x;a)y| < 1 for all yeS. We restrict our attention to 
y_Lx since J(x; a)x = 0. 
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Let y £ E with y_Lx. Then 

y T ((m - l)Ax™- 2 ) y 



|y J J(x;a)y| = 



a 



A + a 



The assumption that (A, x) is negative stable means that C(A, x) is negative definite; 
therefore, y T ((m — l)Ax m ~ 2 ) y < A. On the other hand, by the definition of (3, 

p((m - l)Ax m - 2 ) <(3(A), 

Thus, using the fact that A + a is positive, we have 

Q< -(3(A) + a < y T ((to - l)Ax m - 2 ) y + a A + a ^ 
A + a — A + a A + a 

Hence, p(J(x;a)) < 1, and x is a linearly attracting fixed point. 

On the other hand, if (A,x) is not negative stable, then there exists w £ £ such 
that w_Lx and w T ((to - l)Ax" 1 - 2 ) w > A. Thus, 

T , s w T ((m-l)Jlx m - 2 )w + a X + a 

w T J(x; a)w = ^ 1 1 > = i. 

X + a ~ X + a 

Consequently, p( J(x; a)) > 1, and x is not a linearly attracting fixed point by Theo- 
rem 2.11 and Theorem 2.12. D 

In fact, we can see from the proof of Theorem 4.8 that if the eigenpair (A, x) is 
not negative stable, there is no choice of a £ M. that will make p(J(x.;a)) < 1. For 
x to be a fixed point at all, we must have A + a > 0, and this is sufficient to obtain 
p(J(x; a)) > 1 if (A,x) is not negative stable. In other words, smaller values of a do 
not induce "accidental" convergence to any additional eigenpairs. 

An alternative argument establishes, for a > (3(A), the slightly broader result 
that any attracting fixed point, regardless of order of convergence, must be a strict 
constrained local maximum of /(x) = Ax. m on S. That is, the marginally attract- 
ing case corresponds to a stationary point that has degenerate C(A, x) but is still 
a maximum. This follows from Theorem 2.8, where the needed convexity holds for 
a > (3(A), so that any vector x' £ E in the neighborhood of convergence of x must 
satisfy /(x') < /(x). One can convince oneself that the converse also holds for 
a > (3(A), i.e., any strict local maximum corresponds to an attracting fixed point. 
This is because the strict monotonicity of / under iteration (other than at a fixed 
point) implies that the iteration acts as a contraction on the region of closed contours 
of / around the maximum. 

The counterpart of Theorem 4.8 for the concave case is as follows. 

Corollary 4.9. Let (A,x) be an eigenpair of a symmetric tensor A £ M\ m ' n h 
Assume a £ R such that a < —(3(A), where (3(A) is as defined in (4-2). Let </>(x) be 
given by (4-9). Then (A,x) is positive stable if and only if x is a linearly attracting 
fixed point of —(f). 

Example 4.10. We return again to A £ M^ 3 ! as defined in Example 3.5. 
Figure 4.3a shows the spectral radius of the Jacobian of the fixed point iteration for 
varying values of a for all eigenpairs that are positive or negative stable. At a = 0, 
the spectral radius is greater than 1 for every eigenvalue, and this is why S-HOPM 
never converges. At a — 2, on the other hand, we see that the spectral radius is less 
than 1 for all of the negative stable eigenpairs. Furthermore, the spectral radius stays 
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less than 1 as a increases. Conversely, at a = —2, the spectral radius is less than 1 
for all the eigenpairs that are positive stable. 

In Figure 4.4a, we plot example iteration sequences for ||xfc +1 — x*||/||xfc — x*|| 
for each eigenpair, using a = 2 for the negative stable eigenpairs and a — —2 for 
the positive stable eigenpairs. We expect ||xfc+i — x*|| = <r||xfc — x*|| where a is 
the spectral radius of the Jacobian J(x;a). For example, for A = —1.0954, we have 
a = 0.4 (shown as a dashed line) and this precisely matched the observed rate of 
convergence (shown as a solid line). 

Figure 4.5a plots /(x) on the unit sphere using color to indicate function value. 
We show the front and back of the sphere. Notice that the horizontal axis is from 1 to 
— 1 in the left plot and from —1 to 1 in the right plot, as if walking around the sphere. 
In this image, the horizontal axis corresponds to X2 and the vertical axis to x^; the 
left image is centered at x\ = 1 and the right image at x\ = — 1. Since m is even, 
the function is symmetric, i.e., /(x) = /(— x). The eigenvectors are shown as white, 
gray, and black circles corresponding to their classification as negative stable, positive 
stable, and unstable, respectively; in turn, these correspond to maxima, minima, and 
saddle points of /(x). 

Figure 4.5b shows the basins of attraction for SS-HOPM with a = 2. Every grid 
point on the sphere was used as a starting point for SS-HOPM, and it is colored 4 
according to which eigenvalue it converged to. In this case, every run converges 
to a negative stable eigenpair (labeled with a white circle). Recall that SS-HOPM 
must converge to some eigenpair per Theorem 4.4, and Theorem 4.8 says that it is 
gcnerically a negative stable eigenpair. Thus, the non-attracting points lie on the 
boundaries of the domains of attraction. 

Figure 4.5c shows the basins of attraction for SS-HOPM with a = —2. In this 
case, every starting point converges to an eigenpair that is positive stable (shown as 
gray circles). □ 





-5-4-3-2-1012345 -5-4-3-2-1012345 

a a 

(a) A 6 K [4 ' 3 ) from Example 3.5. (b) A £ R [3 ' 31 from Example 3.6. 

Fig. 4.3: Spectral radii of the Jacobian J(x; a) for different eigenpairs as a varies. 



Example 4.11. We return again to A e Ml 3 ' 3 ! from Example 3.6, which i 



is 



4 Specifically, each block on the sphere is colored according to the convergence of its lower left 
point. 
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(a) A S Ml 4,3 ' from Example 3.5 using a = ±2, (b) .A £ Rl 3 ' 3 ! from Example 3.6 using o = 1. 
as appropriate. 

Fig. 4.4: Example plots of ||xfc+i — x* ||/||xfc — x*||. The expected rate of convergence 
from J(x* ; a) is shown as a dashed line. 



interesting because S-HOPM was able to find 2 of its eigenpairs without any shift. In 
Figure 4.6a, /(x) is plotted on the unit sphere, along with each eigenvector, colored 
white, gray, or black based on whether it is negative stable, positive stable, or unsta- 
ble, respectively. Observe that the function is antisymmetric, i.e., /(x) = — /(— x). 
Figure 4.6b shows the basins of attraction for S-HOPM (i.e., SS-HOPM with a = 0). 
Every starting point converges to one of the 2 labeled eigenpairs. This is not sur- 
prising because Figure 4.3b shows that there are 2 eigenvalues for which the spectral 
radius of the Jacobian is less than 1 (A = 0.8730 and 0.4306). The other 2 eigenvalues 
are non-attracting for a = 0. Figure 4.4b shows the observed rates of convergence. 

Figure 4.6c shows the basins of attraction for SS-HOPM with a = 1; each negative 
stable eigenpair (shown as a white circle) is an attracting eigenpair. The concave case 
is just a mirror image and is not shown. □ 

As the previous example reminds us, for odd order, there is no need to try both 
positive and negative a because the definiteness of C flips for eigenvectors of opposite 
sign. 

Two additional examples of SS-HOPM are presented in Appendix A. 

4.3. Relationship to power method for matrix eigenpairs. The power 
method for matrix eigenpairs is a technique for finding the largest-magnitude eigen- 
value (and corresponding eigenvector) of a diagonalizable symmetric matrix [10]. Let 
A be a symmetric real- valued n x n matrix. Then the matrix power method is defined 
by 



Xfc+l = 



Ax fc 
I Ax* 



Assume that VAV is the Schur decomposition of A with eigenvalues satisfying 
I Ai| > | A 2 1 > •■ > |A„| (note the strict difference in the first 2 eigenvalues). The 
sequence {xfc} produced by the matrix power method always converges (up to sign) 
to the eigenvector associated with Ai- Shifting the matrix by A <— A + al shifts 
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(a) Function values for /(x) = .Ax" 
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(b) SS-HOPM basins of attraction using a = 2. 
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(c) SS-HOPM basins of attraction using a = — 2. 

Fig. 4.5: Illustrations for A G Rl 4 ' 3 ' from Example 3.5. The horizontal axis corre- 
sponds to X2 and the vertical axis to x^; the left image is centered at X\ — 1 and the 
right at X\ = — 1. White, gray, and black dots indicate eigenvectors that are negative 
stable, positive stable, and unstable, respectively. 
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(b) SS-HOPM basins of attraction using a = 0. 
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(c) SS-HOPM basins of attraction using a = 1. 

Fig. 4.6: Illustrations for ./I g Ml 3,3 ! from Example 3.6. The horizontal axis corre- 
sponds to X2 and the vertical axis to X3] the left image is centered at X\ = 1 and the 
right at x\ = — 1. White, gray, and black dots indicate eigenvectors that are negative 
stable, positive stable, and unstable, respectively. 
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the eigenvalues by Xj <— Xj + a, potentially altering which eigenvalue has the largest 
magnitude. 

In the matrix case, the eigenvalues of the Jacobian defined by (4.10) for an eigen- 
pair (AjjXj) are given by 

{0} U I — : 1 < i < n with i ^ j 

[Xj + a 

Thus, the Jacobian at xi is the only one such that p(J(x.;a)) < 1; no other eigen- 
vectors are stable fixed points of the iteration. This corresponds to Theorem 4.8 (or 
Corollary 4.9), since the most positive eigenvalue is negative stable, the most negative 
eigenvalue is positive stable, and every other eigenvalue is unstable. The eigenpair 
(Ai,xi) is an attractor for ordinary (convex) power iteration if Ai > or for flipped 
(concave) power iteration if Ai < 0. 

In contrast to the matrix power method, SS-HOPM can find multiple eigenpairs 
since there may be multiple positive and negative stable eigenpairs. But, as for 
matrices, since the most positive and most negative eigenvalues correspond to the 
global maximum and minimum of /(x), they must be negative stable and positive 
stable respectively. Thus, choosing a positive is necessary for finding the most positive 
tensor eigenvalue; conversely, a negative is necessary for finding the most negative 
tensor eigenvalue. Unfortunately, the ability to find multiple eigenpairs means that 
there is no guarantee that the iterates will converge to an extremal eigenpair from 
every starting point. In fact, multiple starting points may be needed. 

4.4. Comparison to other methods. SS-HOPM is useful for its guaranteed 
convergence properties and its simple implementation based on tensor-vector multi- 
plication. For fixed m and large n, the computational complexity of each iteration 
of SS-HOPM is 0(n m ), which is the number of individual terms to be computed in 
,/tx m_1 . This is analogous to the 0(n 2 ) complexity of matrix-vector multiplication as 
used in the matrix power method. We do not yet know how the number of iterations 
needed for numerical convergence of SS-HOPM depends on m and n. 

The convergence of SS-HOPM to only a subset of eigenvalues, which tend to be 
among the largest in magnitude, is beneficial when the large eigenvalues are of primary 
interest, as in the rank-1 approximation problem [11]. In particular, the most positive 
eigenvalue and most negative eigenvalue always have a region of stable convergence 
for a suitable choice of shift. However, the lack of stable convergence to certain other 
eigenvalues is a disadvantage if those eigenvalues are of interest. 

One evident computational approach for finding tensor eigenpairs should be com- 
pared with SS-HOPM. This is to apply a numerical solver for nonlinear equation 
systems, such as Newton's method, directly to the eigenvalue equations (1.2). The 
computational complexity of each iteration of Newton's method for this system is that 
of SS-HOPM plus the construction and inversion of the (n + 1) X (n + 1) Jacobian 
for (A,x). The Jacobian construction is effectively included in SS-HOPM, since it is 
dominated by computing Ax m ~ 2 , which is a precursor of Ax m ~ l . The additional 
work for inversion is 0(n 3 ), and for m > 3 it does not affect the complexity scaling, 
which remains 0{n m ). 

Two advantages of an approach such as Newton's method are generic locally 
stable convergence, which enables finding eigenpairs not found by SS-HOPM, and the 
quadratic order of convergence, which can be expected to require fewer iterations than 
the linearly convergent SS-HOPM. On the other hand, there is no known guarantee 
of global convergence as there is for SS-HOPM, and it is possible that many starting 
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points fail to converge. Even those that do converge may lead to eigenpairs of less 
interest for a particular application. Furthermore, certain tensor structures can be 
more efficiently handled with SS-HOPM than with Newton's method. For example, 
consider a higher-order symmetric tensor expressed as a sum of terms, each of which is 
an outer product of matrices. The computation of Ax m_1 then reduces to a series of 
matrix-vector multiplications, which are 0(n 2 ). This compares favorably to the 0(n 3 ) 
of Newton's method for the same tensor. Further investigation of general nonlinear 
solver approaches to the tensor eigenvalue problem will be beneficial. 

Finally, we consider a polynomial solver approach, such as we implemented in 
Mathematica. This can find all eigenpairs (subject to numerical conditioning issues) 
but becomes computationally expensive for large m and n. In part this is simply 
because, from Theorem 5.3, the number of eigenpairs grows exponentially with n. 
The solver in Mathematica is designed to find all solutions; it is not clear whether a 
substantial improvement in efficiency would be possible if only one or a few solutions 
were required. 

Nevertheless, for comparison with the iterative approaches discussed above, we 
have measured the computational time per eigenpair on a desktop computer for var- 
ious values of m and n, as shown in Figure 4.7. The complexity of the polynomial 
solution, even measured per eigenpair, is seen to increase extremely rapidly (faster 
than exponentially) with n. Thus the polynomial solver approach is not expected to 
be practical for large n. 



1.00 
0.50 

0.20 

0.10 
0.05 

0.02 

0.01 




Fig. 4.7: Average time (over 10 trials) required to compute all eigenpairs, divided by 
the number of eigenpairs, for random symmetric tensors in ]Ri m >"J. Note logarithmic 
vertical scale. Measured using NSolve in Mathematica on a 4 GHz Intel Core i7. 



5. Complex case. We present the more general definition of complex eigenpairs 
and some related results, and propose an extension of the SS-HOPM algorithm to this 
case. 

5.1. Eigenrings. Some of the solutions of the polynomial system that results 
from the eigenvalue equations may be complex; thus, the definition can be extended 
to the complex case as follows, where t denotes the conjugate transpose. 

Definition 5.1. Assume that A. is a symmetric m -order n- dimensional real- 
valued tensor. Then A G C is an eigenvalue of A. if there exists x £ C" such that 



Ax 7i 



Ax and 



x'x 



1. 



(5.1) 
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The vector*, is a corresponding eigenvector, and (A, x) is called an eigenpair. 

Definition 5.1 is closely related to the E-eigenpairs denned by Qi [18, 19] but 
differs in the constraint on x. 5 It can also be considered as the obvious extension of 
(P-)eigenpairs to C. 

It has been observed [19, 3] that the complex eigenpairs of a tensor form equiva- 
lence classes under a multiplicative transformation. Specifically, if (A, x) is an eigen- 
pair of A £ Rl m,n ] and y = e l¥> x with ip G R, then y f y = x f x = 1 and 

Ay™- 1 = e^^^Ax™ -1 = e' i(m " 1)v Ax = e 4(m ~ 2)v Ay. 

Therefore (e l ' m ~ 2 - )v A, e lv x) is also an eigenpair of A for any <p £ R. Consequently, if 
A is an eigenvalue, then any other A'eC with |A'| = |A| is also an eigenvalue. This 
leads to the notion of an eigenring. 

Definition 5.2 (Eigenring). For any (A, x) £ CxC™ that is an eigenpair 
of A £ Rl m '™] ; we define a corresponding equivalence class of (vector-normalized) 
eigenpairs 

■P(A,x) = {(A',x') : A' = e l(m - 2)v A,x' = e lv x, if £ R}, 

as well as a corresponding eigenring 

K(A) = {A'eC: |A'| = |A|}. 

Thus, any eigenring contains either 1 or 2 real eigenvalues. The special case of real 
eigenpairs occurs whenever the corresponding eigenvector for one of those real eigen- 
values is also real. 

Since we assume that A is real- valued, any nonreal eigenpairs must come in sets 
of 2 related by complex conjugation, because taking the conjugate of the eigenvalue 
equation does not change it. Such conjugate eigenpairs are not members of the same 
equivalence class unless they are equivalent to a real eigenpair. 

An elegant result has recently been derived for the number of distinct (non- 
equivalent) eigenvalues of a symmetric tensor. The result was first derived for even- 
order tensors by Ni et al. [16] and later generalized by Cartwright and Sturmfels [3] 
for all cases. The case of m — 2 requires application of l'Hopital's rule to see that 
there are n eigenvalues. 

Theorem 5.3 (Cartwright and Sturmfels [3]). A generic symmetric tensor A £ 
Rl m > n ] has ((m — l) n — l)/(m — 2) distinct eigenvalue equivalence classes. 

These papers [16, 3] use the condition x T x = 1 to normalize eigenpairs, but in 
the generic case the result is the same for our condition x^x = 1. This is because 
the eigenpairs with x^x = 1 generically consist of isolated equivalence classes that 
have x T x ^ and thus can be rescaled to satisfy x T x = 1, giving a one-to-one 
correspondence between the distinct eigenvalues in the two normalizations. In special 
cases, however, the condition x^x = 1 admits additional eigenpairs with x T x = 0. 
Furthermore, tensors can be constructed with a continuous family of non-equivalent 
eigenvectors that correspond to the same eigenvalue when normalized by x T x but to 
infinitely many distinct eigenvalues when normalized by x^x [3, Example 5.7]. 

The polynomial system solver using Grobner bases mentioned earlier can also be 
used to find complex solutions. A complication is that our normalization condition 
x^x = 1 is nonpolynomial due to the complex conjugation. The system, however, 



5 Qi [18, 19] requires x T x = 1 rather than x^x = 1. 



Shifted Power Method for Computing Tensor Eigenpairs 25 



A 


x' JI 


0.6694 


[ 0.2930 + 0.0571i 0.8171 - 0.0365J -0.4912 - 0.0245J ] 


0.6694 


[ 0.2930 - 0.0571i 0.8171 + 0.0365J -0.4912 + 0.0245i ] 



Tabic 5.1: Nonreal eigenpairs for A <E Rl 4 ' 3 ! from Example 3.5. 



becomes polynomial if the alternate normalization condition x T x = 1 is temporarily 
adopted. Any such x can be rescaled to satisfy x^x = 1. Complex eigenvectors with 
x T x = will not be found, but, as just noted, these do not occur generically. Any 
nonreal solutions are transformed to a representative of the eigenring with positive 
real A by setting 



(A,x) 



|A|\ 1/(m_2) x 



(xtx)™/2-i' \\ J (xtx)V 



This polynomial system solution becomes prohibitively expensive for large m or n; 
however, for Example 3.5, the nonreal eigenpairs can be computed this way and are 
presented in Table 5.1. Thus, from this and Table 3.1, we have found the 13 eigenvalue 
equivalence classes (real and nonreal) guaranteed by Theorem 5.3. 

5.2. SS-HOPM for Complex Eigenpairs. We propose an extension of the 
SS-HOPM algorithm to the case of complex vectors in Algorithm 3. Observe that the 
division by A& + a keeps the phase of x/j from changing unintentionally It is akin to 
taking the negative in the concave case in Algorithm 2. It is important to note that 
even if an eigenpair is real, there is no guarantee that the complex SS-HOPM will 
converge to the real eigenpair; instead, it will converge to some random rotation in 
the complex plane. We have no convergence theory in the convex case, but we present 
several promising numerical examples. 

Algorithm 3 Complex SS-HOPM 

Given a tensor A € R[ m <"] . 

Require: xo £ C n with ||xo|| = 1. Let Ao = Ax™. 
Require: a e C 
1: for k = 0, 1, . . . do 

x fc+ i «- (Ax™- 1 + ax fe )/(Afe + a) 

x fc+ i 4- x fc+ i/||x fc+ i|| 



A fc +i <- xi^Ax™^ 1 



end for 



Example 5.4. We once again revisit A € Rl 4 ' 3 ! from Example 3.5 and test the 
complex version of SS-HOPM in Algorithm 3. Table 5.2a shows the results of 100 
runs using the same experimental conditions as in Example 3.5 except with complex 
random starting vectors. We find 7 distinct eigenrings — the 6 stable real eigenpairs 
as well as a ring corresponding to the 2 complex eigenpairs. Figure 5.1a shows the 
individual A* values plotted on the complex plane. As mentioned above, it may 
converge anywhere on the eigenring, though there is clear bias toward the value of a. 

To investigate this phenomenon further, we do another experiment with a = 
— (1 + i)/y2. It finds the same eigenrings as before as shown in Table 5.2b, but this 
time the A* values are distributed mostly in the lower left quadrant of the complex 
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plane as shown in Figure 5.1b, again close to the value of a. In the case of the 2 
complex eigenpairs with the same eigenring, the method finds the 2 distinct eigenvec- 
tors (i.e., defining 2 different equivalence classes) in the 4 different times it converges 
to that eigenvalue; this is not surprising since the complex eigenvalue has 2 different 
eigenvectors as shown in Table 3.1. 

We also ran an experiment with a — 0. In this case, 95 trials converged, but to 
non-eigenpairs (all with |A| = 0.3656). In each case, even though {A^} converged, we 
had ||xfe+i — Xfc|| — > 1.1993, indicating that the sequence {x^} had not converged and 
hence we did not obtain an x* satisfying the eigenvalue equation (5.1). Although it is 
not shown, in all the convergent examples with the shifts mentioned above, the {x^.} 
sequence converged. □ 



Table 5.2: Eigenrings computed for A € 
HOPM with 100 random starts. 



[[4,3] f rom Example 3.5 by complex SS- 



(a) a = 2 




# Occurrences 


|A| 


18 


1.0954 


18 


0.8893 


21 


0.8169 
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0.6694 


22 


0.5629 


8 


0.3633 


12 


0.0451 



(b) a = v/2(l+i) (2 failures). 
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(b) a = ~\/2(l + i) (2 failures). 



Fig. 5.1: For A. S M} 4 ' 3 ^ from Example 3.5, final A values (indicated by red asterisks) 
for 100 runs of complex SS-HOPM. The green lines denote the eigenrings. 



6. Conclusions. We have developed the shifted symmetric higher-order power 
method (SS-HOPM) for finding tensor eigenvalues. The method can be considered 
as a higher-order analogue to the power method for matrices. Just as in the matrix 
case, it cannot find all possible eigenvalues, but it is guaranteed to be able to find the 
largest-magnitude eigenvalue. Unlike the matrix case, it can find multiple eigenvalues; 
multiple starting points are typically needed to find the largest eigenvalue. A GPU- 
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based implementation of SS-HOPM has been reported [1]. 

Building on previous work [11, 23, 8], we show that SS-HOPM will always converge 
to a real eigenpair for an appropriate choice of a. Moreover, using fixed point analysis, 
we characterize exactly which real eigenpairs can be found by the method, i.e., those 
that are positive or negative stable. Alternative methods will need to be developed 
for finding the unstable real eigenpairs, i.e., eigenpairs for which C(A,x) is indefinite. 
A topic for future investigation is that the boundaries of the basins of attraction for 
SS-HOPM seem to be defined by the non-attracting eigenvectors. 

We present a complex version of SS-HOPM and limited experimental results indi- 
cating that it finds eigenpairs, including complex eigenpairs. Analysis of the complex 
version is a topic for future study. 

Much is still unknown about tensor eigenpairs. For example, how do the eigen- 
pairs change with small perturbations of the tensor entries? Is there an eigendecom- 
position of a tensor? Can the convergence rate of the current method be accelerated? 
How does one numerically compute unstable eigenpairs? For computing efficiency, 
what is the optimal storage for symmetric tensors? These are all potential topics of 
future research. 
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Appendix A. Further examples. For additional insight, we consider two 
analytical examples. In the experiments in this section, each random trial used a point 
Xo chosen from a uniform distribution on [—1,1]". We allow up to 1000 iterations 
and say that the algorithm has converged if |Afc+i — Afc| < 10 -15 . 

Consider the tensor A. € K^ 3 ' 3 ' whose entries are except where the indices are 
all unequal, in which case the entries are 1, i.e., 



Qijk 



1 if (i, j, k) is some permutation of (1,2,3), 
otherwise. 



(A.l) 



Any eigenpair (A, x) must satisfy the following system of equations: 



2x 2 x 2 = Xxi, 



2xix 3 = Xx 2 , 



2x\x 2 — AX3 



xj + xi + xj = 1. 



The 7 real eigenpairs can be computed analytically and are listed in Table A. la, from 
which we can see that there are 4 negative stable eigenpairs that should be identifiable 
by SS-HOPM. Figure A. la shows the spectral radius of the Jacobian as a varies; the 
curve is identical for all 4 negative stable eigenpairs. 

Another example is given as follows. Define the tensor A. <E R[ 4,2 1 by 

dijkl — for all i,j,k,l except ami = 1 and a 2222 = —I. (A. 2) 

The eigenpairs can be computed analytically as solutions to the following system: 



x\ 



\xi, 



Xx-2 



x\ 



xl = \. 



The 2 real eigenpairs are listed in Table A. lb, from which we can see that one is 
negative stable and the other is positive stable. Figure A. la shows the spectral radius 
of the Jacobian as a varies. In this case, the spectral radius of the Jacobian can be 
computed analytically; for A — 1, it is jx— and hence there is a singularity for a = — 1. 



Table A.l: Eigenpairs for two analytical problems. 

(a) Eigenpairs for A £ IRl 3 ' 3 ' from (A.l). 



A 


X 


eigenvalues of C(A, x) 


Type 





[100] 


{-2,2} 


Unstable 





[010] 


{-2,2} 


Unstable 





[001] 


{-2,2} 


Unstable 


2/V3 


[ i/Va i/%/3 


1A/3] 


{ -2.3094 , -2.3094 } 


Neg. Stable 


2/V3 


[ l/Vs-i/Vs- 


-l/v/3] 


{ -2.3094 , -2.3094 } 


Neg. Stable 


2/V3 


[-1A/3 I/n/3- 


-l/v/3] 


{ -2.3094 , -2.3094 } 


Neg. Stable 


2/V3 


[-1A/3-1/V3 


1/V3] 


{ -2.3094 , -2.3094 } 


Neg. Stable 



(b) Eigenpairs for A e M [4 ' 21 from (A.2). 



A 


X 


eigenvalues of C(A, x) 


Type 


1 


[10] 


{-1} 


Neg. Stable 


-1 


[01] 


{1} 


Pos. Stable 



For (A.l), we ran 100 trials with a = 0, and none converged, as expected per 
Figure A. la. The results of 100 random trials with a = 12 (the "conservative choice") 
is shown in Table A. 2a, in which case every trial converged to one of the 4 negative 
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(a) (A.l) (All four are identical.) (b) (A.2) 

Fig. A.l: Spectral radii of the Jacobian J(x.;a) for different eigenpairs as a varies. 
Table A.2: Eigenpairs for A G R [3 > 3] from (A.l) computed by SS-HOPM. 







(a) a = 12 






# Occurrences 


A 


X 


Median Its. 


22 


1.1547 


[ -0.5774 0.5774 


-0.5774 ] 


92 


18 


1.1547 


0.5774 0.5774 


0.5774 ] 


90 


29 


1.1547 


-0.5774 -0.5774 


0.5774 ] 


91 


31 


1.1547 


0.5774 -0.5774 


-0.5774 ] 


94 







(b) a = 1 






# Occurrences 


A 


X 


Median Its. 


22 


1.1547 


[ 0.5774 -0.5774 


-0.5774 ] 


9 


25 


1.1547 


[ -0.5774 0.5774 


-0.5774 ] 


9 


26 


1.1547 


0.5774 0.5774 


0.5774 ] 


9 


27 


1.1547 


-0.5774 -0.5774 


0.5774 ] 


9 



stable eigenpairs. (Note that 2/y/3 « 1.1547 and l/y/3 « 0.5774.) Table A. 2b shows 
the results of 100 random trials with a = 1. As expected (per Figure A. la), the 
convergence is much faster. For (A.2), we ran 100 trials with a = 0.5 (Table A. 3a) 
and 100 trials with a — —0.5 (Tabic A. 3b). We find the negative stable and positive 
stable eigenvalues as expected. 



Table A.3: Eigenpairs for A € R [4 ' 2] from (A.2) computed by SS-HOPM. 





(a) a = 0.5 




# Occurrences 


A 


X 


Median Its. 


100 


1.0000 


[ -1.0000 0.0000 ] 


16 





(b)a = 


= -0.5 






# Occurrences 


A 


X 


Median Its. 


100 


-1.0000 


[ 


- 0.0000 


1.0000 ] 


16 



